---
  title: "CARRS RDD analysis :Robustness check keeping only follow uo observations within 12 months"
knit: (function(input, ...) {rmarkdown::render(input, output_file="")})
output: pdf_document
geometry: "left=2cm,right=2cm,top=1cm,bottom=1cm"
date: "`r format(Sys.time(), '%d %B, %Y')`"
---

```{r , setup, include=FALSE, eval = TRUE}
knitr::opts_chunk$set(
   comment = '', warning=FALSE,  message=FALSE, results= FALSE , fig.width = 9,fig.height = 7, echo=FALSE, root.dir =""

)

 knitr::opts_knit$set(
   root.dir = ""    
 )
```

```{r packages}
library(tidyverse)
library(gridExtra)
library(dplyr)
library(rdrobust)
library(rddensity)
library(viridisLite)
library(foreign)
library(ggplot2)
library(data.table)
library(doBy)
library(socviz)
library(openxlsx) 
library(modelsummary)
library(labelled)
library(gtsummary)
library(gt)
library(kableExtra)
```

```{r loaddata, eval = TRUE}
# data for outcome: drug treatment
  data_diagnosis<-data.frame(read.csv("carrs_selected.csv"))
  data_diagnosis<-as.data.table(data_diagnosis)
 
# data for outcome: diagnosis
  data_treatment<-data.frame(read.csv("carrs_full.csv"))
  data_treatment<-as.data.table(data_treatment)
  
# data for outcomes: systolic and diastolic blood pressure
  CARRSbp <-as.data.table(read.csv("carrs_bp.csv"))
  CARRSbp <- na.omit(CARRSbp, cols=(c("sysdiff","diasdiff","edu")))
```

## binary outcomes
##calculating the bandwidth for each sub group

``` {r RunningVariableAll, eval = TRUE}
rv.update = function(b,c) {
  
  out=c[city==b,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                            (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
  
  
  a_b<<-out
  
}

rv.update("Delhi",data_diagnosis)
diag_delhi<-a_b[city=="Delhi",]
rm(a_b)

rv.update("Chennai",data_diagnosis)
diag_chennai<-a_b[city=="Chennai",]
rm(a_b)


rv.update("Delhi",data_treatment)
treat_delhi<-a_b[city=="Delhi",]
rm(a_b)

rv.update("Chennai",data_treatment)
treat_chennai<-a_b[city=="Chennai",]
rm(a_b)

```

``` {r Running VariableBySex, eval = TRUE}
# by sex
rv.update = function(a,b,c) {
  
  out=c[female==a & city==b,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                            (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
  
  
  a_b<<-out
  
}

# males
rv.update(0,"Delhi",data_diagnosis)
diag_delhiM<-a_b[city=="Delhi" & female==0,]
rm(a_b)

rv.update(0,"Chennai",data_diagnosis)
diag_chennaiM<-a_b[city=="Chennai" & female==0,]
rm(a_b)


rv.update(0,"Delhi",data_treatment)
treat_delhiM<-a_b[city=="Delhi" & female==0,]
rm(a_b)

rv.update(0,"Chennai",data_treatment)
treat_chennaiM<-a_b[city=="Chennai" & female==0,]
rm(a_b)

# females
rv.update(1,"Delhi",data_diagnosis)
diag_delhiF<-a_b[city=="Delhi" & female==1,]
rm(a_b)

rv.update(1,"Chennai",data_diagnosis)
diag_chennaiF<-a_b[city=="Chennai" & female==1,]
rm(a_b)


rv.update(1,"Delhi",data_treatment)
treat_delhiF<-a_b[city=="Delhi" & female==1,]
rm(a_b)

rv.update(1,"Chennai",data_treatment)
treat_chennaiF<-a_b[city=="Chennai" & female==1,]
rm(a_b)
```

# RDD function - binary outcomes

```{r RDfunction, eval = TRUE }
##Defining the common function##
carrs_rdd = function(a,b,d) {
  
  d$Y<-a
  d$X<-b
  
  
  out =  rdrobust(y = d$Y, d$X, kernel = "triangular", p = 1, covs = cbind(d$w01, d$w23, d$w45),
                  c = 0, cluster = d$pid, rho = 1 , all = T)
  summary(out)
  
  ret <- data.frame(
    City = c(formatC(round(100*out$coef[1,1],2),2,format="f"), formatC(round(out$pv[3,1],2),2,format="f"), 
            paste0( "(", formatC(round(100*out$ci[3,1],2),2,format="f"), " - ", formatC(round(100*out$ci[3,2],2),2,format="f"), ")" ), 
            formatC(round(out$bws[1,1],2),2,format="f"), as.character(out$N_h[1] + out$N_h[2]))
  )
  rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  
  output<<-ret %>% add_rownames  
}
```

# Outcome: diagnosed in the past 12 months
```{r Diagnosis, eval = TRUE}
rdd = function(a) {
  
  carrs_rdd(a$hypdiag,a$runningvar,a)
  
}

rdd(diag_delhi)
diag_delhi <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(diag_chennai)
diag_chennai <-output%>%mutate(Chennai=City)%>%select(-City)

# males
rdd(diag_delhiM)
diag_delhiM <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(diag_chennaiM)
diag_chennaiM <-output%>%mutate(Chennai=City)%>%select(-City)

# females
rdd(diag_delhiF)
diag_delhiF <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(diag_chennaiF)
diag_chennaiF <-output%>%mutate(Chennai=City)%>%select(-City)
```

# Outcome: current treatment
```{r Treatment, eval = TRUE}
rdd = function(a) {
  
  carrs_rdd(a$hypdrug,a$runningvar,a)
  
}

rdd(treat_delhi)
treat_delhi <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(treat_chennai)
treat_chennai <-output%>%mutate(Chennai=City)%>%select(-City)

# males
rdd(treat_delhiM)
treat_delhiM <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(treat_chennaiM)
treat_chennaiM <-output%>%mutate(Chennai=City)%>%select(-City)

# females
rdd(treat_delhiF)
treat_delhiF <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(treat_chennaiF)
treat_chennaiF <-output%>%mutate(Chennai=City)%>%select(-City)
```


```{r MergingResults, eval = TRUE}
diag_table<-merge(diag_chennai,diag_delhi,by='rowname')
diag_table<-diag_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

treat_table<-merge(treat_chennai,treat_delhi,by='rowname')
treat_table<-treat_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

# males
diag_tableM<-merge(diag_chennaiM,diag_delhiM,by='rowname')
diag_tableM<-diag_tableM %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

treat_tableM<-merge(treat_chennaiM,treat_delhiM,by='rowname')
treat_tableM<-treat_tableM %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

# females
diag_tableF<-merge(diag_chennaiF,diag_delhiF,by='rowname')
diag_tableF<-diag_tableF %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

treat_tableF<-merge(treat_chennaiF,treat_delhiF,by='rowname')
treat_tableF<-treat_tableF %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))
```

```{r Table, eval = TRUE}
  filename <- paste0("city_combined.tex")
  caption=paste0("Effect of home-based screening on hypertension diagnosis and treatment by city")

  comb_diag_table<-rbind(diag_table,diag_tableF,diag_tableM)
  comb_treat_table<-rbind(treat_table,treat_tableF,treat_tableM)
  
  comb_prevdiag_table <- cbind(comb_diag_table,comb_treat_table)
  comb_prevdiag_table <- comb_prevdiag_table[-c(4)]

combined<-kbl(comb_prevdiag_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccc", format = "latex",col.names = c(" ","Chennai","Delhi","Chennai","Delhi"),digits = 2) %>%
    add_header_above(c(" ", "Diagnosis" = 2, "Treatment" = 2)) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Females", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10) %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Males", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15) %>%
    footnote(general = c("Note: The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and confidence intervals resulting from robust bias-corrected standard errors clustered at the individual-level.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T) %>%
save_kable(filename, format = "latex")
```

## continuous outcomes
## calculating the bandwidth for each subgroup - continuous outcomes

``` {r RunningVariableAll, eval = TRUE}
rv.update = function(b,c) {
  
  out=c[city==b,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                            (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
   a_b<<-out
  
}


rv.update("Delhi",CARRSbp)
bp_delhi<-a_b[city=="Delhi",]
rm(a_b)

rv.update("Chennai",CARRSbp)
bp_chennai<-a_b[city=="Chennai",]
rm(a_b)
```

``` {r Running VariableBySex, eval = TRUE}
# by sex
rv.update = function(a,b,c) {
  
  out=c[female==a & city==b,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                            (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
  
  
  a_b<<-out
}

# males
rv.update(0,"Delhi",CARRSbp)
bp_delhiM<-a_b[city=="Delhi" & female==0,]
rm(a_b)

rv.update(0,"Chennai",CARRSbp)
bp_chennaiM<-a_b[city=="Chennai" & female==0,]
rm(a_b)

# females
rv.update(1,"Delhi",CARRSbp)
bp_delhiF<-a_b[city=="Delhi" & female==1,]
rm(a_b)

rv.update(1,"Chennai",CARRSbp)
bp_chennaiF<-a_b[city=="Chennai" & female==1,]
rm(a_b)
```


# RDD function - continuous outcomes

```{r RDfunctionBP, eval = TRUE }
##Defining the common function##
carrs_rdd = function(a,b,d) {
  
  d$Y<-a
  d$X<-b
  
  
  out =  rdrobust(y = d$Y, d$X, kernel = "triangular", p = 1, covs = cbind(d$w02, d$w24),
                  c = 0, cluster = d$pid, rho = 1 , all = T)
  summary(out)
  
  ret <- data.frame(
    City = c(formatC(round(out$coef[1,1],2),2,format="f"), formatC(round(out$pv[3,1],2),2,format="f"), 
            paste0( "(", formatC(round(out$ci[3,1],2),2,format="f"), " - ", formatC(round(out$ci[3,2],2),2,format="f"), ")" ), 
            formatC(round(out$bws[1,1],2),2,format="f"), formatC(out$N_h[1] + out$N_h[2]))
  )
  rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  
  output<<-ret %>% add_rownames  
}
```

# Outcome: systolic
```{r Systolic, eval = TRUE}
rdd = function(a) {
  
  carrs_rdd(a$sysdiff,a$runningvar,a)
  
}

rdd(bp_delhi)
sys_delhi <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(bp_chennai)
sys_chennai <-output%>%mutate(Chennai=City)%>%select(-City)

# males
rdd(bp_delhiM)
sys_delhiM <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(bp_chennaiM)
sys_chennaiM <-output%>%mutate(Chennai=City)%>%select(-City)

# females
rdd(bp_delhiF)
sys_delhiF <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(bp_chennaiF)
sys_chennaiF <-output%>%mutate(Chennai=City)%>%select(-City)
```

# Outcome: current diastolic
```{r Diastolic, eval = TRUE}
rdd = function(a) {
  
  carrs_rdd(a$diasdiff,a$runningvar,a)
  
}

rdd(bp_delhi)
dias_delhi <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(bp_chennai)
dias_chennai <-output%>%mutate(Chennai=City)%>%select(-City)

# males

rdd(bp_delhiM)
dias_delhiM <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(bp_chennaiM)
dias_chennaiM <-output%>%mutate(Chennai=City)%>%select(-City)

# females

rdd(bp_delhiF)
dias_delhiF <-output%>%mutate(Delhi=City)%>%select(-City)

rdd(bp_chennaiF)
dias_chennaiF <-output%>%mutate(Chennai=City)%>%select(-City)
```


```{r MergingResults, eval = TRUE}
sys_table<-merge(sys_chennai,sys_delhi,by='rowname')
sys_table<-sys_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

dias_table<-merge(dias_chennai,dias_delhi,by='rowname')
dias_table<-dias_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

# males
sys_tableM<-merge(sys_chennaiM,sys_delhiM,by='rowname')
sys_tableM<-sys_tableM %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

dias_tableM<-merge(dias_chennaiM,dias_delhiM,by='rowname')
dias_tableM<-dias_tableM %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

# females
sys_tableF<-merge(sys_chennaiF,sys_delhiF,by='rowname')
sys_tableF<-sys_tableF %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

dias_tableF<-merge(dias_chennaiF,dias_delhiF,by='rowname')
dias_tableF<-dias_tableF %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))
```

```{r Table, eval = TRUE}
  filename <- paste0("city_combined_bpdiff.tex")
  caption=paste0("Effect of home-based screening on systolic and diastolic blood pressure by city")


  comb_sys_table_bp<-rbind(sys_table,sys_tableF,sys_tableM)
  comb_dias_table_bp<-rbind(dias_table,dias_tableF,dias_tableM)
  
  comb_citybp_table <- cbind(comb_sys_table_bp,comb_dias_table_bp)
  comb_citybp_table <- comb_citybp_table[-c(4)]

combined<-kbl(comb_citybp_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccc", format = "latex",col.names = c(" ","Chennai","Delhi","Chennai","Delhi"),digits = 2)  %>%
    add_header_above(c(" "=1, "Systolic" = 2, "Diastolic" = 2)) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Females", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10) %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Males", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15) %>%
    footnote(general = c("Note: The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and confidence intervals resulting from robust bias-corrected standard errors clustered at the individual-level.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T) %>%
save_kable(filename, format = "latex")
```
